1 Análise de séries temporais

No curso pretendemos explorar vários conceitos teóricos relacionados à séries temporais, passando pela definição formal, propriedades estatísticas, distribuições finito dimensionais, modelos, entre outros conceitos. No entanto, pela minha experiência trabalhando na prática, fazendo acessorias e também desenvolvendo metodologias para análise de séries temporais, uma deficiência que notei foi a dificuldade de lidar com características dos dados que não são explicadas por componentes aleatórias, e sim por componentes determinísticas dos dados, como a tendência e a sazonalidade. Esse material é dedicado a explorar a ler uma série temporal no R, identificar diferentes tipos de tendências e sazonalidades e apresentar diferentes formas de lidar com essas questões de uma forma prática, sem maiores preocupações coma formalidade conceitual.

A maior parte da teoria estatística clássica é baseada na suposição de aleatoriedade da amostra e observações independentes. Por outro lado, na série temporal é exatamente o oposto. As observações geralmente são dependentes de valores anteriores, e sua análise deve levar em conta sua ordem temporal.

Por exemplo, um estudo de prospective cohort comparando taxa de lesões antes e depois de um programa implementado, a análise de séries temporais considera a variabilidade que ocorre ao longo do período de estudo além da mudança associada à intervenção . Isso evita a perda de informações sobre a variabilidade da incidência ao longo do tempo, quando as taxas são agregadas em uma taxa antes e outra depois. A população é seu próprio controle.

Se observações anteriores podem prever observações futuras com exatidão, temos um processo determinístico. No entanto, a previsão exata geralmente é impossível, pois os valores passados determinam apenas parte do valor futuro. Neste caso, dizemos que o processo é estocástico, e os valores futuros devem ser vistos como uma probabilidade condicionada aos valores passados.

Existem muitos modelos disponíveis para descrever o comportamento de uma determinada série. A escolha de tal modelo depende de fatores como o comportamento do fenômeno ou o conhecimento prévio de sua natureza e a finalidade da análise.

1.1 O que é uma série temporarl?

Dados de séries temporais são obeservações de um evento ou fenômeno ao longo do tempo. Os intervalos de observações pode ou não ser igualmente espaçados. Geralmente são anos, trimestres, meses, semanas, dias, horas, minutos e segundos. Mas outros tipos de espaçamentos entre observações também são comuns. Como é o caso do Censo Demográfico.

Na Figura 1.1 apresentamos alguns exemplos de séries temporais com diferentes intervalos de observações.

Figura 1.1: Séries temporais com diferentes intervalos de observações

Nesse exemplo temos séries temporais com diferentes intervalos de tempo. A série histórica do Pib tem intervalos de tempo de um ano. Os dados sobre valor da cesta básica em Porto alegre é divulgado mensalmente, os dados sobre Covid19 são divulgados diariamente e uma série temporal com observações a cada dez anos é a da população brasileira.

library(MASS)
data() #rodar no seu próprio R-studio

1.1.1 Diferentes tipos de séries temporais

A instalação R base já nos fornece muitos conjuntos de dados para trabalhar em séries temporais. Para este exemplo, primeiro carregamos o pacote MASS que contém alguns dos conjuntos de dados que usaremos. Podemos listar todos os conjuntos de dados disponíveis com a função data() do pacote utils.

1.1.1.1 Séries temporais medidas em intervalos de tempo regulares (discretos).

Econômico: bolsa de valores

plot(EuStockMarkets[, 1], main = "Preços de fechamento diários do DAX", ylab = "Preço")

Físico/Biológico: Pluviometria, DNA

plot(((nhtemp) - 32) * 5 / 9, main = "Temperaturas médias anuais em New Haven", ylab = "Temperatura em ºC")

Marketing: vendas por mês

plot(BJsales, main = "Dados de vendas diárias (150 dias)", ylab = "Vendas")

Demografia: população por ano, acidentes de carro por dia

plot(austres, main = "Série temporal trimestral do número de residentes australianos", ylab = "Residentes")

Controle de processo: medições de fábrica como pesos finais das latas, índices de qualidade

library(qicharts2)

## Warning: package 'qicharts2' was built under R version 4.1.3

qic(((nhtemp) - 32) * 5 / 9,  chart = 'i')

1.1.1.2 Medidas em intervalos de tempo irregulares (eventos)

Processos pontuais: terremotos (eventos)

eq_data <- readr::read_csv("https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_day.csv")
eq_data$yes <- 0
plot(eq_data$time, eq_data$yes, main = "Terremotos de magnitude 2,5+ USGS, no dia anterior",
     ylab = "Evento", xlab = "Hora do dia", yaxt = 'n', sub = paste("Dados de ", Sys.Date()))

1.1.1.3 Medidas continuamente

Processos binários: teoria da comunicação (liga/desliga, zeros e uns);

set.seed(2114)
binary_data <- unlist(replicate(10, rep(runif(1) < 0.5, floor(runif(1, 10, 20)))))
binary_data <- stepfun(seq_len(length(binary_data) - 1), binary_data)
plot(binary_data, main = "Exemplo binário sintético", ylim = c(0,1),
     ylab = "Value", xlab = "Tempo contínuo", do.points = FALSE, yaxp = c(0, 1, 1))

Sinais analógicos: som, temperatura, umidade, ECG

ecg<- read.csv("https://raw.githubusercontent.com/marciovalk/timeseries_booklet/main/data/ECG.csv",skip=1)
plot.ts(ecg,xlab="",ylab="tempo")

1.2 Quais os objetivos da análise de séries temporais?

O estudo de séries temporais tem por objetivos principais definir o processo gerador de dados, fazer previsões futuras da série, identificar ciclos, tendências e/ou sazonalidades de forma que a decisão que envolve as variáveis em questão seja a mais acurada possível.

Dada uma série temporal $\{Z(t_1 ),\ldots, Z(t_N)\}$, observada nos instantes $t_1,\ldots,t_N$, podemos estar interessados em:

1. Investigar o mecanismo gerador da série temporal;
1. Fazer previsões de valores futuros da série; podendo ser a curto ou longo prazo;
1. Descrever apenas o comportamento da série através de gráficos;
1. Procurar periodicidades relevantes nos dados. Em todos estes casos podemos construir modelos probabilísticos ou estocásticos, tanto no domínio do tempo como no domínio da freqüência, por exemplo: um sinal aleatório com frequência medida em Hz. Devemos construir modelos simples e com menor número de parâmetros possíveis.

Uma série temporal pode ser decomposta em seus componentes de forma a entendê-la, analisá-la, modelá-la e fazer previsões sistematicamente. Esta é uma introdução para iniciantes à análise de séries temporais, respondendo a perguntas fundamentais como:

Quais são os componentes de uma série temporal?
O que é uma série temporal estacionária?
Como fazer uma decomposição?
Como lidar com a tendência?
Como dessazonalizar uma série temporal?
O que é correlação? etc.

Séries temporais serão abordadas ao longo do texto, no qual não exploramos as propriedades teóricas das séries, mas sim as características amostrais, sempre com um enfoque prático usando o R.

1.3 Como criar uma série temporal no R?

O inputData é idealmente um vetor numérico da classe numeric ou integer. Depois de ler os dados de série temporal no R, a próxima etapa é armazenar os dados em um objeto da classe ts, para que você possa usar as muitas funções do R para analisar dados de séries temporais. Para armazenar os dados em um objeto de série temporal, usamos a função ts(). Por exemplo, para armazenar os dados na variável ‘kings’ como um objeto de série temporal em R, digitamos:

#ts(inputData, start=c(2009), end=c(2014), frequency=1) # Dados anuais
#ts(inputData, frequency = 4, start = c(1959, 2)) # frequency 4 => Dados trimestrais
time_data1 <- ts(1:10, frequency = 12, start = 1990) # freq 12 => Dados mensais. 
time_data1

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 1990   1   2   3   4   5   6   7   8   9  10

class(time_data1)

## [1] "ts"

1.4 Lendo uma série temporal no R

Você pode ler dados em R usando a função scan(), que assume que seus dados para pontos de tempo sucessivos estão em um arquivo de texto simples com uma coluna.

Por exemplo, o arquivo https://github.com/marciovalk/timeseries_booklet/blob/main/data/kings_age.txt contém dados sobre a idade de morte de sucessivos reis da Inglaterra, começando com William, o Conquistador (fonte original: Hipel e Mcleod, 1994). O conjunto de dados fica assim:

Figura: Idade da morte de sucessivos reis da Inglaterra

Apenas as primeiras linhas do arquivo foram mostradas. As duas primeiras linhas contêm algum comentário sobre os dados, e queremos ignorar isso ao lermos os dados em R. Podemos fazer usando o parâmetro skip da função scan(), que especifica quantas linhas no topo do arquivo devem ser ignoradas. Ler o arquivo em R, ignorando as duas primeiras linhas, pode ser feito da seguinte forma:

kings_data <- scan("https://raw.githubusercontent.com/marciovalk/timeseries_booklet/main/data/kings_age.txt",skip=2)

Neste caso, a idade da morte de 42 reis sucessivos da Inglaterra foi lida na variável kings_data. A próxima etapa é armazenar os dados em um objeto do tipo ts em R, para que você possa usar as muitas funções para analisar dados de séries temporais. Podemos usar a função ts() para armazenar os dados na variável kings_data como um objeto do tipo ts em R:

 kingstimeseries <- ts(kings_data)
 kingstimeseries

## Time Series:
## Start = 1 
## End = 42 
## Frequency = 1 
##  [1] 60 43 67 50 56 42 50 65 68 43 65 34 47 34 49 41 13 35 53 56 16 43 69 59 48
## [26] 59 86 55 68 51 33 49 67 77 81 67 71 81 68 70 77 56

Às vezes, o conjunto de dados de série temporal pode ter sido coletado em intervalos regulares inferiores a um ano, por exemplo, mensalmente ou trimestralmente. Nesse caso, podemos especificar o número de vezes que os dados foram coletados por ano usando o parâmetro frequency na função ts(). Para dados de séries temporais mensais, define-se frequency=12, enquanto para dados de séries temporais trimestrais, o input frequency=4.

O primeiro ano em que os dados foram coletados e o primeiro intervalo nesse ano deve ser indicado no parâmetro start na função ts(). Por exemplo, se o primeiro ponto de dados corresponde ao segundo trimestre de 1986, define-se start=c(1986,2).

Um exemplo é um conjunto de dados do número de nascimentos por mês na cidade de Nova York, de janeiro de 1946 a dezembro de 1959 (originalmente coletados por Newton). Esses dados estão disponíveis no arquivo https://github.com/marciovalk/timeseries_booklet/blob/main/data/nybirths.csv Podemos ler os dados no R e armazená-los como um objeto ts, digitando:

nyb<- read.csv("https://raw.githubusercontent.com/marciovalk/timeseries_booklet/main/data/nybirths.csv")
nybtimeseries <- ts(nyb, frequency=12, start=c(1946,1))
nybtimeseries

##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 1946 23.598 26.931 24.740 25.806 24.364 24.477 23.901 23.175 23.227 21.672
## 1947 21.089 23.709 21.669 21.752 20.761 23.479 23.824 23.105 23.110 21.759
## 1948 20.035 23.590 21.672 22.222 22.123 23.950 23.504 22.238 23.142 21.059
## 1949 20.000 22.424 20.615 21.761 22.874 24.104 23.748 23.262 22.907 21.519
## 1950 20.894 24.677 23.673 25.320 23.583 24.671 24.454 24.122 24.252 22.084
## 1951 23.049 25.076 24.037 24.430 24.667 26.451 25.618 25.014 25.110 22.964
## 1952 22.270 24.775 22.646 23.988 24.737 26.276 25.816 25.210 25.199 23.162
## 1953 22.644 25.565 24.062 25.431 24.635 27.009 26.606 26.268 26.462 25.246
## 1954 23.304 26.982 26.199 27.210 26.122 26.706 26.878 26.152 26.379 24.712
## 1955 24.239 26.721 23.475 24.767 26.219 28.361 28.599 27.914 27.784 25.693
## 1956 24.218 27.914 26.975 28.527 27.139 28.982 28.169 28.056 29.136 26.291
## 1957 24.848 27.543 26.896 28.878 27.390 28.065 28.141 29.048 28.484 26.634
## 1958 24.924 28.963 26.589 27.931 28.009 29.229 28.759 28.405 27.945 25.912
## 1959 25.286 27.660 25.951 26.398 25.565 28.865 30.000 29.261 29.012 26.992
##         Nov    Dec
## 1946 21.870 21.439
## 1947 22.073 21.937
## 1948 21.573 21.548
## 1949 22.025 22.604
## 1950 22.991 23.287
## 1951 23.981 23.798
## 1952 24.707 24.364
## 1953 25.180 24.657
## 1954 25.688 24.990
## 1955 26.881 26.217
## 1956 26.987 26.589
## 1957 27.735 27.132
## 1958 26.619 26.076
## 1959 27.897

1.5 Plotando séries temporais

Depois de ler uma série temporal no R, o próximo passo geralmente é fazer um gráfico dos dados da série temporal, o que você pode fazer com a função plot.ts().

Por exemplo, para traçar a série temporal da idade da morte de 42 reis sucessivos da Inglaterra, podemos fazer da seguinte forma:

plot.ts(kingstimeseries)

Podemos ver no gráfico que esta série temporal provavelmente poderia ser descrita usando um modelo aditivo, em que cada ponto de dados ($Y_t$) no tempo $t$ em uma Série Temporal pode ser expresso como uma soma de 3 componentes

\[Y_t=S_t+T_t+A_t\] em que $S_t$ é a sazonalidade, $T_t$ é a tendência e $A_t$ é uma componente aleatória. No modelo aditivo, as flutuações aleatórias nos dados são aproximadamente constantes em tamanho ao longo do tempo.

Da mesma forma, para traçar a série temporal do número de nascimentos por mês na cidade de Nova York, digitamos:

plot.ts(nybtimeseries)

Podemos ver nesta série temporal que parece haver variação sazonal no número de nascimentos por mês: há um pico a cada verão e um vale a cada inverno. Novamente, parece que esta série temporal provavelmente pode ser descrita usando um modelo aditivo, pois as flutuações sazonais são aproximadamente constantes em tamanho ao longo do tempo e não parecem depender do nível da série temporal, e as flutuações aleatórias também parecem ser aproximadamente constante em tamanho ao longo do tempo.

No entanto, na série temporal das vendas mensais da loja de souvenirs em uma cidade em Queensland, Austrália, representada na seguinte Figura,

sou <- scan("https://raw.githubusercontent.com/marciovalk/timeseries_booklet/main/data/souvenir.txt")
souts <- ts(sou, frequency=12, start=c(1987,1))
plot.ts(souts)

o modelo aditivo não não aparenta ser adequado para descrever esta série temporal, uma vez que o tamanho das flutuações sazonais e as flutuações aleatórias parecem aumentar com o nível da série temporal. Assim, podemos precisar transformar a série temporal para obter uma série temporal transformada que possa ser descrita usando um modelo aditivo. Por exemplo, podemos transformar a série temporal utilizando o logaritmo natural dos dados originais:

logsouts<-log(souts)
plot.ts(logsouts)

1.6 Outros exemplos de séries temporais

Nem sempra as características dos dados são observáveis no gráfico da série temporal. Por exemlo, o Southern Oscillation Index (SOI) é um índice padronizado baseado nas diferenças observadas de pressão ao nível do mar (SLP) entre Tahiti e Darwin, Austrália. Os dados são encontrado no site https://www.ncei.noaa.gov/access/monitoring/enso/soi, ou podem ser obtidos diretamente do pacote ocedata.

if(!require("ocedata")) install.packages("ocedata")

## Warning: package 'ocedata' was built under R version 4.1.3

require("ocedata")
data(soi, package="ocedata")

plot(soi$year, soi$index, type='l', xlab="Year", ylab="SOI")

Podemos observar os dados anuais, flutuando em torno de uma média constante. Além disso, não é muito claro qualquer outra característica dos dados.

Na série temporal de Temperatura Global, em que dados são incluídos a partir da análise de Surface Temperature (GISTEMP) e do componente global do Climate at a Glance (GCAG), dois conjuntos de dados são fornecidos: 1) média mensal global e 2) anomalias de temperatura média anual em graus Celsius de 1880 até o presente. https://datahub.io/core/global-temp#readme Nesse caso os dados são estruturados de uma forma que facilita a utilização desses no formato data.frame, em que é possível utilizar as funções do pacote ggplot2.

if(!require("plotly")) install.packages("plotly")
require(plotly)
require(ggplot2)
gtemp<- read.csv('https://datahub.io/core/global-temp/r/annual.csv')

p <- ggplot(data=gtemp, aes(x=Year, y=Mean, group=Source)) +
    geom_line(color="#00BFC4") +
    geom_point(color="#F8766D")

fig <- ggplotly(p)

fig

Nese exemplo, temos uma característica de tendência não linear, em que a média dos dados aparenta ter um crescimento quadrático.

2 Antes da modelagem de séries temporais

Para a parte inferencial da análise de séries temporais algumas suposições, comuns em várias metodologias estatísticas, também são feitas nesse contexto. Entre elas estão a homogeneidade de variância e a normalidade dos dados. No entanto, essas propriedades nem sempre são satisfeitas pelos dados, sendo necessárias algumas transformações.

2.1 Transformações para estabilizar a variância

Para estabilizar a variância, uma transformação logarítmica pode ser a solução. Outra possibilidade é uma raíz quadrada dos dados.

opar <- par(mfrow = c(2, 1), mar = c(3.1, 4.1, 2.1, 1.1))

claims <- ts(Insurance$Claims, start = c(1973, 365.25*3/4), frequency = 365.25)
plot.ts(claims, main = "Números de sinistros de seguro de carro (1973)", ylab = "Reivindicações")
par(mar = c(5.1, 4.1, 0.1, 1.1))
plot(log(claims), ylab = "log(Reivindicações)")

par(opar)

Uma transformação logarítmica transforma uma sazonalidade multiplicativa em aditiva. No entanto, isso só estabilizará a variância se o termo de erro também for multiplicativo.

2.2 Transformações para normalizar os dados

Os dados são geralmente considerados normais. Transformação logarítmica e raiz quadrada podem ser usadas, no entanto, são apenas casos especiais da transformação Box-Cox.

library(graphics)
layout(mat = matrix(c(1, 1, 2, 3), ncol = 2, byrow = TRUE)) # setting plot 
plot.ts(airquality$Ozone, main = "Medições de qualidade do ar de Nova York, maio a setembro de 1973", ylab = "Ozônio", xlab = "Tempo")

ozone <- ts(na.omit(airquality$Ozone))
hist(ozone, probability = 1); lines(density(ozone), col = "red") # distribuição dos dados
hist(log(ozone), probability = 1); lines(density(log(ozone)), col = "red") #  normalizado?

É interessante notar que Nelson e Granger, em 1979, encontraram pouco benefício na aplicação de uma transformação geral Box-Cox em vários conjuntos de dados. Normalmente, é aconselhável aplicar o mínimo de transformações possíveis, exceto quando a variável tem uma interpretação física direta.

2.3 Decomposição de séries temporais

Uma forma de começar a analisar uma série temporal é separando-a em seus componentes, o que significa separá-la em componente de tendência, uma componente sazonal/cíclica e uma componente aleatória.

2.3.1 Decompondo dados com tendência- abordagem paramétrica

Podemos ajustar uma reta através da função lm() do R em então retirar a tendência fazendo a série real menos a reta ajustada.

oldpar <- par(mfrow = c(2, 1), mar = c(3.1, 4.1, 2.1, 1.1))

plot.ts(kingstimeseries)

fit <- lm((kingstimeseries) ~ seq_along(kingstimeseries))
pred <- predict(fit, data.frame(seq_along(kingstimeseries)))
pred <- ts(pred, start = start(kingstimeseries), frequency = frequency(kingstimeseries))
# `pred` contains our line based on the simple linear model we did.
data <- kingstimeseries-pred # why divide and not subtract? because the variance also increases (we have a multiplicative seasonality)
lines(pred, col = "red")
par(mar = c(5.1, 4.1, 0.1, 1.1))
plot(data, ylab = "kings age - fitted")
par(oldpar)

2.3.2 Abordagens não-paramétrica para lidar com a tendência

Também chamadas de Filtering, essas abordagens são operações um pouco mais complexa do que o ajuste de curvas. Não estamos apenas tentando encontrar um polinômio que melhor se ajuste aos nossos dados, mas estamos transformando nossa TS original em outra TS usando uma fórmula (que aqui chamamos de filtro). Este filtro pode ser um dos vários tipos de Médias móveis, regressões ponderadas localmente (por exemplo, LOESS) ou Splines (um polinômio por partes). Uma ressalva dessas técnicas de suavização é o problema do efeito final (já que em uma extremidade da série temporal não temos todos os valores para calcular, por exemplo, a média móvel)

\[SMA(x_t)=\frac{1}{2q+1}\sum_{r=-q}^{+q}x_{t+r}.\]

2.3.2.1 Média móvel simples

Uma série temporal não sazonal consiste em uma componente de tendência e uma componente aleatória. De forma bem geral, poderíamos escrever da forma

\[Y_t=\mu_t+\varepsilon_t\] em que $mu_t$ descreve a média(tendência) ao longo do tempo, e $\varepsilon_t$ é uma componente aleatória/estocástica com alguma distribuição de probabilidade.

Decompor a série temporal envolve tentar separar a série temporal nesses componentes, ou seja, estimar o componente de tendência e o componente irregular.

Para estimar o componente de tendência de uma série temporal não sazonal que pode ser descrita usando um modelo aditivo, é comum usar algum método de suavização, como a média móvel simples da série temporal.

A função SMA() no pacote R TTR pode ser usada para suavizar dados de séries temporais usando uma média móvel simples.

if(!require("TTR")) install.packages("TTR")
library("TTR")

A função SMA() pode ser usada para suavizar uma série temporal, basta especificar a ordem (span) da média móvel simples, usando o parâmetro “n”. Por exemplo, para calcular uma média móvel simples de ordem 5, definimos n=5 na função SMA(). Por exemplo, como discutido acima, a série temporal da idade da morte de 42 reis sucessivos da Inglaterra parece não ser sazonal

kingstimeseriesSMA5 <- SMA(kingstimeseries,n=5)
plot.ts(kingstimeseries)
plot.ts(kingstimeseriesSMA5)

Ainda parece haver muitas flutuações aleatórias na série temporal suavizada usando uma média móvel simples de ordem 5. Assim, para estimar o componente de tendência com mais precisão, podemos tentar suavizar os dados com uma média móvel simples de um ordem superior.

 kingstimeseriesSMA15 <- SMA(kingstimeseries,n=15)
plot.ts(kingstimeseries)
plot.ts(kingstimeseriesSMA15)

No entanto, observe que estamos “perdendo” as n primeiras informações. Encontrar o n agequado para suavização é uma tarefa de tentativa e erro e depende muito do contexto e da experiência do pesquisador.

2.3.2.2 Splines

Outra maneira de “encaixar uma linha” é usando um spline. O número de nós é sua escolha, mas quanto menor, mais suave (próximo á uma reta) é a curva.

oldpar <- par(mfrow = c(3, 1), mar = c(3.1, 4.1, 2.1, 1.1))
plot(USAccDeaths, main = "Mortes Acidentais entre 1973-1978 nos EUA", ylab = "Mortes Acidentais")

pred <- smooth.spline(USAccDeaths, nknots = 5)$y
pred <- ts(pred, start = start(USAccDeaths), frequency = frequency(USAccDeaths))
lines(pred, col = "red")
par(mar = c(5.1, 4.1, 0.1, 1.1))
data <- USAccDeaths - pred 
plot(data, ylab = "Deaths - spline")
par(mar = c(5.1, 4.1, 0.1, 1.1))
data <- USAccDeaths / pred
plot(data, ylab = "Deaths / spline")

par(oldpar)

2.3.3 Diferenciação

É um tipo especial de filtragem, onde calculamos a diferença entre o valor atual e o próximo (retorno). É útil para remover tendências, tornando uma TS estacionária (sem tendência). Podemos usar a diferenciação várias vezes (chamamos ordem), mas geralmente uma iteração (primeira ordem) é suficiente. O operador matemático usado para denotar a diferenciação é o nabla ($\nabla$). $\nabla^2$ significa diferenciação de segunda ordem.

oldpar <- par(mfrow = c(2, 1), mar = c(3.1, 4.1, 2.1, 1.1))
plot(co2, main = "Concentração atmosférica de CO2 de Mauna Loa", ylab = "concentração de CO2 ")
par(mar = c(5.1, 4.1, 0.1, 1.1))
plot(diff(co2), ylab = "diff(co2)")

par(oldpar)

A primeira diferença remove toda a tendência, mas ainda temos a sazonalidade.

Cuidado (Efeito Slutzky-Yule) Eles mostraram que, usando operações como diferenciação e média móvel, pode-se induzir variação senoidal nos dados que, na verdade, não são informações reais.

2.4 Decompondo dados sazonais

Uma série temporal sazonal consiste de uma componente de tendência, uma componente sazonal e uma componente aleatória. Decompor a série temporal significa separar a série temporal nesses três componentes: ou seja, estimar esses três componentes.

Para estimar o componente de tendência e o componente sazonal de uma série temporal sazonal que pode ser descrita usando um modelo aditivo, podemos usar a função decompose() em R. Esta função estima os componentes tendência, sazonal e aleatória de uma série temporal que pode ser descrito usando um modelo aditivo.

A função decompose() retorna um objeto de lista como resultado, onde as estimativas do componente sazonal, componente de tendência e componente aleatória são armazenadas em elementos nomeados desses objetos de lista, chamados de seasonal, trend, and random respectivamente.

Por exemplo, como discutido acima, a série temporal do número de nascimentos por mês na cidade de Nova York é sazonal com pico a cada verão e vale a cada inverno, e provavelmente pode ser descrita usando um modelo aditivo, uma vez que as flutuações sazonais e aleatórias parecem ser aproximadamente constante em tamanho ao longo do tempo:

A tendência, componente sazonal e componente aleatória pode ser estimada da forma:

nybts_comp <- decompose(nybtimeseries)

Os valores estimados dos componentes sazonal, de tendência e aleatória agora são armazenados nas variáveis nybts_comp$seasonal*, *nybts_comp$trend e nybts_comp$random. Por exemplo, podemos imprimir os valores estimados do componente sazonal digitando:

nybts_comp$seasonal

##             Jan        Feb        Mar        Apr        May        Jun
## 1946 -2.0893513  0.8561327 -0.8080692  0.2452609 -0.1596461  1.5263413
## 1947 -2.0893513  0.8561327 -0.8080692  0.2452609 -0.1596461  1.5263413
## 1948 -2.0893513  0.8561327 -0.8080692  0.2452609 -0.1596461  1.5263413
## 1949 -2.0893513  0.8561327 -0.8080692  0.2452609 -0.1596461  1.5263413
## 1950 -2.0893513  0.8561327 -0.8080692  0.2452609 -0.1596461  1.5263413
## 1951 -2.0893513  0.8561327 -0.8080692  0.2452609 -0.1596461  1.5263413
## 1952 -2.0893513  0.8561327 -0.8080692  0.2452609 -0.1596461  1.5263413
## 1953 -2.0893513  0.8561327 -0.8080692  0.2452609 -0.1596461  1.5263413
## 1954 -2.0893513  0.8561327 -0.8080692  0.2452609 -0.1596461  1.5263413
## 1955 -2.0893513  0.8561327 -0.8080692  0.2452609 -0.1596461  1.5263413
## 1956 -2.0893513  0.8561327 -0.8080692  0.2452609 -0.1596461  1.5263413
## 1957 -2.0893513  0.8561327 -0.8080692  0.2452609 -0.1596461  1.5263413
## 1958 -2.0893513  0.8561327 -0.8080692  0.2452609 -0.1596461  1.5263413
## 1959 -2.0893513  0.8561327 -0.8080692  0.2452609 -0.1596461  1.5263413
##             Jul        Aug        Sep        Oct        Nov        Dec
## 1946  1.1582032  0.6852257  0.7688539 -1.1161557 -0.3832102 -0.6835852
## 1947  1.1582032  0.6852257  0.7688539 -1.1161557 -0.3832102 -0.6835852
## 1948  1.1582032  0.6852257  0.7688539 -1.1161557 -0.3832102 -0.6835852
## 1949  1.1582032  0.6852257  0.7688539 -1.1161557 -0.3832102 -0.6835852
## 1950  1.1582032  0.6852257  0.7688539 -1.1161557 -0.3832102 -0.6835852
## 1951  1.1582032  0.6852257  0.7688539 -1.1161557 -0.3832102 -0.6835852
## 1952  1.1582032  0.6852257  0.7688539 -1.1161557 -0.3832102 -0.6835852
## 1953  1.1582032  0.6852257  0.7688539 -1.1161557 -0.3832102 -0.6835852
## 1954  1.1582032  0.6852257  0.7688539 -1.1161557 -0.3832102 -0.6835852
## 1955  1.1582032  0.6852257  0.7688539 -1.1161557 -0.3832102 -0.6835852
## 1956  1.1582032  0.6852257  0.7688539 -1.1161557 -0.3832102 -0.6835852
## 1957  1.1582032  0.6852257  0.7688539 -1.1161557 -0.3832102 -0.6835852
## 1958  1.1582032  0.6852257  0.7688539 -1.1161557 -0.3832102 -0.6835852
## 1959  1.1582032  0.6852257  0.7688539 -1.1161557 -0.3832102

Os fatores sazonais estimados são dados para os meses de janeiro a dezembro e são os mesmos para cada ano. O maior fator sazonal é para julho (cerca de 1,46), e o menor é para fevereiro (cerca de -2,08), indicando que parece haver um pico nos nascimentos em julho e um mínimo nos nascimentos em fevereiro de cada ano.

plot(nybts_comp$seasonal)

O gráfico acima mostra a série temporal original (topo), o componente de tendência estimada (segundo a partir do topo), a componente sazonal estimada (terceiro a partir do topo) e a componente aleatória estimada (abaixo). Vemos que a componente de tendência estimada mostra um pequeno decréscimo de cerca de 24 em 1947 para cerca de 22 em 1948, seguido por um aumento constante a partir de então para cerca de 27 em 1959.

2.4.1 Ajuste sazonal

Se você tiver uma série temporal sazonal que possa ser descrita usando um modelo aditivo, poderá ajustar sazonalmente a série temporal estimando o componente sazonal e subtraindo o componente sazonal estimado da série temporal original. Podemos fazer isso usando a estimativa do componente sazonal calculado pela função decompose().

Por exemplo, para ajustar sazonalmente a série temporal do número de nascimentos por mês na cidade de Nova York, podemos estimar o componente sazonal usando decompose() e, em seguida, subtrair o componente sazonal da série temporal original

nybtimeseriesseasonallyadjusted <- nybtimeseries - nybts_comp$seasonal
plot.ts(nybtimeseries) 
plot(nybtimeseriesseasonallyadjusted)

Você pode ver que a variação sazonal foi removida da série temporal com ajuste sazonal. A série temporal ajustada sazonalmente agora contém apenas o componente de tendência e um componente irregular.

2.4.2 Lidando com a sazonalidade usando diferenças

A sazonalidade também pode ser eliminada diferenciando com defasagem. Por exemplo, com dados mensais, podemos usar o operador $\nabla_{12}$:

\[\nabla_{12}x_t=x_t−x_{t−12}\]

oldpar <- par(mfrow = c(2, 1), mar = c(3.1, 4.1, 2.1, 1.1))
plot(nybtimeseries)
par(mar = c(5.1, 4.1, 0.1, 1.1))

plot(diff(nybtimeseries, lag = 12), ylab = "diff(lag = 12)")

par(oldpar)

Veja que o exemplo acima removeu a sazonalidade, mas manteve a tendência.

3 Previsões usando suavização exponencial (Exponential Smoothing)

A suavização exponencial pode ser usada para fazer previsões de curto prazo para dados de séries temporais.

Aqui usamoso método conhecido com Holt-Winters, cuja função do R é HoltWinters() e os detalhes podem ser acessados no link https://stat.ethz.ch/R-manual/R-patched/library/stats/html/HoltWinters.html .

A suavização exponencial de Holt-Winters estima o nível, a inclinação e a componente sazonal no ponto de tempo atual. A suavização é controlada por três parâmetros: alfa, beta e gama, para as estimativas do nível, inclinação b da componente de tendência e a componente sazonal, respectivamente, no ponto de tempo atual. Os parâmetros alfa, beta e gama têm valores entre 0 e 1, e valores próximos de 0 significam que um peso relativamente pequeno é dado às observações mais recentes ao fazer previsões de valores futuros.

3.1 Suavização exponencial de Holt-Winters (para dados com tendência)

Em séries temporais que podem ser descrita usando um modelo aditivo com tendência crescente ou decrescente e sem sazonalidade, você pode usar a suavização exponencial de Holt-Winters para fazer previsões de curto prazo.

Um exemplo de série temporal que provavelmente pode ser descrita usando um modelo aditivo com tendência e sem sazonalidade é a série temporal do diâmetro anual da bainha das saias femininas, de 1866 a 1911. Os dados estão disponíveis no arquivo https://raw.githubusercontent.com/marciovalk/timeseries_booklet/main/data/skirts.csv (dados originais de Hipel e McLeod, 1994).

skirts <- read.csv("https://raw.githubusercontent.com/marciovalk/timeseries_booklet/main/data/skirts.csv",skip=1,sep=",")
skirtsseries <- ts(skirts$X608,start=c(1866))
plot.ts(skirtsseries)

Podemos ver pelo gráfico que houve um aumento no diâmetro da bainha de cerca de 600 em 1866 para cerca de 1050 em 1880, e que depois o diâmetro da bainha diminuiu para cerca de 520 em 1911.

Para fazer previsões, podemos ajustar um modelo preditivo usando a função HoltWinters() definindo o parâmetro gamma=FALSE (o parâmetro gamma é usado para suavização exponencial de Holt-Winters , como descrito abaixo).

skirtsseriesforecasts <- HoltWinters(skirtsseries, gamma=FALSE)
skirtsseriesforecasts

## Holt-Winters exponential smoothing with trend and without seasonal component.
## 
## Call:
## HoltWinters(x = skirtsseries, gamma = FALSE)
## 
## Smoothing parameters:
##  alpha: 0.8382467
##  beta : 1
##  gamma: FALSE
## 
## Coefficients:
##         [,1]
## a 529.307759
## b   5.689884

O valor estimado de alfa é 0.84 e de beta é 1.00. Ambos são altos, nos dizendo que tanto a estimativa do valor atual do nível quanto a inclinação b do componente de tendência se baseiam principalmente em observações muito recentes na série temporal. Isso faz sentido intuitivo, já que o nível e a inclinação da série temporal mudam bastante ao longo do tempo.

Podemos plotar a série temporal original como uma linha preta, com os valores previstos como uma linha vermelha em cima dela, digitando:

plot(skirtsseriesforecasts)

Podemos ver na imagem que as previsões dentro da amostra concordam muito bem com os valores observados, embora tendam a ficar um pouco atrás dos valores observados.

Se desejar, você pode especificar os valores iniciais do nível e a inclinação b do componente de tendência usando os argumentos l.start e b.start para a função HoltWinters(). É comum definir o valor inicial do nível para o primeiro valor da série temporal (608 para os dados) e o valor inicial da inclinação para o segundo valor menos o primeiro valor (9 para os dados).

HoltWinters(skirtsseries, gamma=FALSE, l.start=608, b.start=9)

## Holt-Winters exponential smoothing with trend and without seasonal component.
## 
## Call:
## HoltWinters(x = skirtsseries, gamma = FALSE, l.start = 608, b.start = 9)
## 
## Smoothing parameters:
##  alpha: 0.8295494
##  beta : 1
##  gamma: FALSE
## 
## Coefficients:
##         [,1]
## a 529.236680
## b   5.644025

3.1.1 Previsões fora da amostra

Para fazer previsões para tempos futuros não incluídos na série temporal original, usamos a função forecast() no pacote forecast.

if(!require("forecast")) install.packages("forecast")
library("forecast")
library("stats")
pred_skirt <- forecast(skirtsseries,h= 8)
plot(pred_skirt)

As previsões são mostradas como uma linha azul, com os intervalos de previsão de 80% e 95% como áreas sombreadas em cinza.

3.2 Suavização Exponencial Holt-Winters (para dados sazonais e sem tendência)

Uma série temporal que provavelmente pode ser descrita usando um modelo aditivo com sazonalidade e sem tendência é a série temporal a temperatura diária da cidade de Porto Alegre-RS O arquivo https://github.com/marciovalk/timeseries_booklet/blob/main/data/TempPortoAlegreA.xls contém a temperatura diária da cidade de Porto Alegre-RS, a partir de 02/01/2011 até 23/06/2020.

if(!require("httr")) install.packages("httr")
library(httr)

github_link <- "https://github.com/marciovalk/timeseries_booklet/raw/main/data/TempPortoAlegreA.xls"

temp_file <- tempfile(fileext = ".xls")
req <- GET(github_link, 
          # authenticate using GITHUB_PAT
           authenticate(Sys.getenv("GITHUB_PAT"), ""),
          # write result to disk
           write_disk(path = temp_file))
data_temp<-readxl::read_excel(temp_file)
poa_temp<-ts(data_temp$TEMPERATURA,frequency = 365,start = c(2011,02,01))
plot.ts(poa_temp)

Para fazer previsões, podemos ajustar um modelo preditivo usando a função HoltWinters(). Por exemplo, para ajustar um modelo preditivo para a temperatura diária da cidade de Porto Alegre-RS, procedemos da seguinte forma:

poa_temp_forecasts <- HoltWinters(poa_temp,beta=FALSE)
poa_temp_forecasts

## Holt-Winters exponential smoothing without trend and with additive seasonal component.
## 
## Call:
## HoltWinters(x = poa_temp, beta = FALSE)
## 
## Smoothing parameters:
##  alpha: 0.8166299
##  beta : FALSE
##  gamma: 1
## 
## Coefficients:
##               [,1]
## a     26.787537430
## s1    -5.245162488
## s2    -6.314524494
## s3    -3.381265585
## s4     0.209120096
## s5     1.964818608
## s6    -0.091847748
## s7    -1.746645070
## s8    -3.790162536
## s9    -6.870137700
## s10   -9.397227001
## s11  -10.209122887
## s12  -12.543219791
## s13  -12.366377177
## s14  -10.236495485
## s15   -7.281155108
## s16   -4.185702121
## s17   -2.779666545
## s18   -2.510844696
## s19   -3.164754618
## s20   -2.360512361
## s21   -1.567221158
## s22   -3.721337896
## s23   -4.841861855
## s24   -3.678915389
## s25   -4.800306245
## s26   -5.550104458
## s27   -7.626992038
## s28   -6.215475524
## s29   -4.163897726
## s30   -2.984213155
## s31   -2.447252292
## s32   -2.807764342
## s33   -4.752332063
## s34   -3.851081148
## s35   -4.104257940
## s36   -4.358912018
## s37   -5.339529486
## s38   -5.685495201
## s39   -8.827361981
## s40   -9.294417139
## s41   -9.985925539
## s42   -7.594887309
## s43   -4.652209655
## s44   -2.926304517
## s45   -3.702945435
## s46   -3.877752638
## s47   -3.104799453
## s48   -0.822053975
## s49   -0.294038367
## s50   -2.052132215
## s51   -1.528168887
## s52    0.681152805
## s53    1.118897257
## s54    0.828640513
## s55   -3.929473881
## s56   -7.064708735
## s57   -9.463255572
## s58   -9.942325135
## s59   -8.197549512
## s60   -5.650725391
## s61   -4.362781692
## s62   -5.376314632
## s63   -5.060678990
## s64   -4.192031954
## s65   -3.710349011
## s66   -1.782175848
## s67   -4.604840049
## s68   -8.442972624
## s69   -7.318013476
## s70   -5.228686749
## s71   -3.767876299
## s72   -3.771933730
## s73   -4.554472560
## s74   -3.777855526
## s75   -2.970200366
## s76   -2.069656452
## s77   -2.719656838
## s78   -4.398627810
## s79   -5.479802736
## s80   -4.655327124
## s81   -3.695793474
## s82   -2.948539546
## s83   -2.014248115
## s84   -1.933232529
## s85   -1.794864951
## s86   -1.640649255
## s87   -1.941269318
## s88   -3.485760931
## s89   -2.978071244
## s90   -2.314938612
## s91   -3.902787480
## s92   -2.713012494
## s93   -2.706555981
## s94   -1.045287804
## s95   -0.273603385
## s96    1.100373767
## s97    1.547788675
## s98    0.975722272
## s99    1.360947662
## s100  -2.918122459
## s101  -3.178977729
## s102  -0.811199050
## s103   0.421553057
## s104   0.398545033
## s105   2.856597262
## s106   3.237768672
## s107   2.159568380
## s108   0.555740330
## s109   0.759443788
## s110   1.638828143
## s111   0.682738853
## s112  -0.701576994
## s113   0.331125363
## s114   0.927491451
## s115   0.541717662
## s116   0.083754200
## s117  -0.722476133
## s118  -0.311694787
## s119   1.488180101
## s120   3.188285401
## s121   2.833574975
## s122   1.974051740
## s123   1.348248161
## s124   0.762246779
## s125  -0.709096762
## s126  -0.873572242
## s127  -0.926032368
## s128  -0.902791612
## s129  -1.425311812
## s130  -1.835328218
## s131   0.277609412
## s132  -0.003446124
## s133   0.690819170
## s134   2.798390218
## s135   2.801075250
## s136   4.169482351
## s137   4.980430477
## s138   4.311962056
## s139   2.112187423
## s140   0.788905971
## s141   0.081574046
## s142   0.477226033
## s143  -0.247717274
## s144   0.576668448
## s145   0.488989577
## s146   1.316652077
## s147   2.327537765
## s148   2.716097550
## s149   3.342198658
## s150   2.918634633
## s151   1.587602505
## s152   1.561037237
## s153   1.894236502
## s154   3.692672335
## s155   5.514229202
## s156   6.427079396
## s157   5.487867791
## s158   4.988222829
## s159   3.320496165
## s160   1.196160774
## s161   1.041031559
## s162   2.374394558
## s163   4.130430002
## s164   3.853561221
## s165   4.791578265
## s166   4.885926670
## s167   3.907182184
## s168   3.894559026
## s169   4.596585388
## s170   4.051056589
## s171   2.918851327
## s172   1.798475012
## s173   2.196884603
## s174   2.212578031
## s175   2.566047806
## s176   3.902489887
## s177   3.210060983
## s178   5.481342241
## s179   6.230203215
## s180   6.192155693
## s181   5.189298654
## s182   3.826449887
## s183   2.952598721
## s184   3.245725058
## s185   3.566891761
## s186   2.894422790
## s187   3.085845595
## s188   3.128386489
## s189   2.807193495
## s190   1.853387560
## s191   1.556300090
## s192   1.197890961
## s193   1.667644499
## s194   2.585352630
## s195   4.732178009
## s196   5.699356586
## s197   4.812507938
## s198   5.384852996
## s199   7.067576463
## s200   7.060150785
## s201   5.200705828
## s202   4.307787482
## s203   3.388688057
## s204   3.034350116
## s205   3.703149564
## s206   2.829759630
## s207   3.394096079
## s208   5.301420374
## s209   5.588532671
## s210   6.325218927
## s211   6.743705235
## s212   8.351317377
## s213   9.240984209
## s214   7.118662535
## s215   4.084078605
## s216   3.721945538
## s217   3.381744767
## s218   4.626341095
## s219   5.271086464
## s220   5.878058304
## s221   5.943411360
## s222   6.489067397
## s223   7.524442479
## s224   6.955261366
## s225   6.876639683
## s226   7.050628540
## s227   7.096212003
## s228   7.007739570
## s229   6.729260004
## s230   5.678424599
## s231   4.194551652
## s232   4.122342584
## s233   4.168903828
## s234   5.187851754
## s235   7.433084616
## s236   8.999235542
## s237  10.213054087
## s238  10.229250269
## s239  10.324140279
## s240   8.955953471
## s241   7.305788305
## s242   5.697598405
## s243   4.895635833
## s244   5.595767059
## s245   5.306664212
## s246   3.843825442
## s247   2.708989282
## s248   3.946283979
## s249   5.233059110
## s250   5.951240289
## s251   5.789982838
## s252   4.437261043
## s253   4.629317035
## s254   5.275990029
## s255   5.588072484
## s256   6.607113035
## s257   6.392788422
## s258   6.206084003
## s259   6.612066110
## s260   6.844914829
## s261   6.804887395
## s262   6.662309338
## s263   5.727506850
## s264   4.254388799
## s265   2.487536782
## s266   1.021644822
## s267   3.127175767
## s268   3.142570749
## s269   2.759372261
## s270   1.817917455
## s271   1.905450769
## s272   2.694484897
## s273   2.483814334
## s274   0.909643678
## s275  -0.127534952
## s276  -0.693645032
## s277  -0.800500166
## s278  -1.034269161
## s279  -0.792521904
## s280  -0.042727078
## s281   1.100543664
## s282   2.088690024
## s283   3.196592763
## s284   1.999134905
## s285   1.558077813
## s286   0.761089998
## s287  -0.162283805
## s288  -0.103137255
## s289   0.512728628
## s290   3.109941826
## s291   2.847156682
## s292   1.773251028
## s293   0.775933647
## s294   1.163085485
## s295   0.470948151
## s296   0.246211739
## s297   1.152350004
## s298   1.189712172
## s299   2.174884790
## s300   2.359358257
## s301  -0.081165004
## s302  -2.379374484
## s303  -3.050468972
## s304  -3.867060351
## s305  -4.254591984
## s306  -4.670880204
## s307  -5.796470599
## s308  -5.841826920
## s309  -5.584376882
## s310  -6.596138789
## s311  -6.910717337
## s312  -5.164043419
## s313  -2.728788803
## s314  -1.751616368
## s315  -0.636077415
## s316  -0.585888942
## s317  -0.666812739
## s318   0.164367845
## s319   1.688368033
## s320   1.396595278
## s321  -0.304403432
## s322  -2.280376435
## s323  -5.012261389
## s324  -4.884585068
## s325  -5.483824966
## s326  -5.052332487
## s327  -4.002175852
## s328  -2.170332968
## s329  -1.553048788
## s330  -1.091666978
## s331  -0.936741858
## s332  -1.380689861
## s333  -2.348638607
## s334  -2.063661155
## s335  -0.982363049
## s336  -0.097946514
## s337   0.042240195
## s338   0.584748682
## s339   1.419793381
## s340  -0.437887400
## s341  -3.225194354
## s342  -4.984498707
## s343  -6.111686088
## s344  -6.116475023
## s345  -6.063250724
## s346  -6.808495842
## s347  -8.953319408
## s348 -11.639788641
## s349 -12.967903102
## s350 -12.860410824
## s351 -10.397287708
## s352  -7.776818641
## s353  -5.255327562
## s354  -2.791706150
## s355  -2.413513489
## s356  -2.675714879
## s357  -2.046763771
## s358  -2.854856655
## s359  -4.476392779
## s360  -6.971662183
## s361  -8.081809465
## s362  -4.997529572
## s363  -5.121970114
## s364  -6.160867938
## s365  -5.477237430

Os valores estimados de alfa e gamma são 0.81 e 1, respectivamente. O valor de alfa (0.81) é relativamente alto, indicando que a estimativa do nível no momento atual é baseada em fortemente em observações recentes e em algumas observações no passado mais distante. O valor de gamma (0.96) é alto, indicando que a estimativa do componente sazonal no momento atual é baseada apenas em observações muito recentes.

Quanto à suavização exponencial simples e à suavização exponencial de Holt, podemos plotar a série temporal original como uma linha preta, com os valores previstos como uma linha vermelha em cima disso:

plot(poa_temp_forecasts)

Vemos no gráfico que o método exponencial de Holt-Winters é muito bem-sucedido em prever os picos sazonais, que ocorrem aproximadamente em janeiro de cada ano.

3.2.1 Previsões fora da amostra

Por exemplo, os dados originais terminam em 23/06/2020. Se quiséssemos fazer previsões para mais um ano, poderíamos fazer da seguinte forma:

pred_temp1 <- forecast(poa_temp_forecasts,h= 365)
plot(pred_temp1)

pred_temp2 <- predict(poa_temp_forecasts, 365, prediction.interval = TRUE,level = 0.90)
plot(poa_temp_forecasts,pred_temp2)

Podemos investigar se o modelo preditivo pode ser melhorado verificando se os erros de previsão na amostra mostram autocorrelações diferentes de zero nas defasagens 1-20, fazendo um correlograma e realizando o teste Ljung-Box:

acf(pred_temp1$residuals, lag.max=20,na.action = na.omit)

Box.test(pred_temp1$residuals, lag=20, type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  pred_temp1$residuals
## X-squared = 755.78, df = 20, p-value < 2.2e-16

O correlograma mostra que as autocorrelações para os erros de previsão na amostra excedem os limites de significância para defasagens 1-20. Além disso, o p-valor para o teste Ljung-Box é bem pequeno indicando que há muita evidência de autocorrelações diferentes de zero nas defasagens 1-20.

Outro passo para verificar se o modelo preditivo pode ser melhorado é investigar se os erros de previsão são normalmente distribuídos com média zero e variância constante. Portanto, faça um gráfico de tempo dos erros de previsão na amostra para visualizar a variância e trace um histograma dos erros de previsão (com uma curva normal sobreposta que tem média zero e o mesmo desvio padrão que a distribuição dos erros de previsão) para verificar para distribuição normal.

plot.ts(pred_temp1$residuals)

##define an R function “plotForecastErrors()”
plotForecastErrors <- function(forecasterrors)
  {
    # make a histogram of the forecast errors:
  mybinsize <- IQR(forecasterrors, na.rm = TRUE)/4
    mysd   <- sd(forecasterrors, na.rm = TRUE)
    mymin  <- min(forecasterrors, na.rm = TRUE) - mysd*5
    mymax  <- max(forecasterrors, na.rm = TRUE) + mysd*3
    # generate normally distributed data with mean 0 and standard deviation mysd
    mynorm <- rnorm(10000, mean=0, sd=mysd)
    mymin2 <- min(mynorm, na.rm = TRUE)
    mymax2 <- max(mynorm, na.rm = TRUE)
    if (mymin2 < mymin ) { mymin <- mymin2}
    if (mymax2 > mymax) { mymax <- mymax2}
    # make a red histogram of the forecast errors, with the normally distributed data overlaid:
    mybins <- seq(mymin, mymax, mybinsize)
    hist(forecasterrors, col="red", freq=FALSE, breaks=mybins)
    # freq=FALSE ensures the area under the histogram = 1
    # generate normally distributed data with mean 0 and standard deviation mysd
    myhist <- hist(mynorm, plot=FALSE, breaks=mybins)
    # plot the normal curve as a blue line on top of the histogram of forecast errors:
    points(myhist$mids, myhist$density, type="l", col="blue", lwd=2)
}
#plot a histogram (with overlaid normal curve) of the forecast errors for the rainfall predictions
plotForecastErrors(pred_temp1$residuals)

O teste Ljung-Box revela forte evidência de autocorrelações diferentes de zero nos erros de previsão na amostra, as variações nos erros de previsão parecem aproximadamente constantes e os erros de previsão parecem distribuídos normalmente. É necessário uma modelagem mais complexa (paramétrica?) para capturar todas as informações que estão nos dados.

3.3 Suavização Exponencial Holt-Winters (para dados sazonais, com tendência)

Um exemplo de uma série temporal que provavelmente pode ser descrita usando um modelo aditivo com tendência e sazonalidade é a série temporal do log das vendas mensais da loja de souvenirs em uma cidade balneária em Queensland, Austrália (logsouts)

souforecasts <- HoltWinters(logsouts)
souforecasts

## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = logsouts)
## 
## Smoothing parameters:
##  alpha: 0.413418
##  beta : 0
##  gamma: 0.9561275
## 
## Coefficients:
##            [,1]
## a   10.37661961
## b    0.02996319
## s1  -0.80952063
## s2  -0.60576477
## s3   0.01103238
## s4  -0.24160551
## s5  -0.35933517
## s6  -0.18076683
## s7   0.07788605
## s8   0.10147055
## s9   0.09649353
## s10  0.05197826
## s11  0.41793637
## s12  1.18088423

Os valores estimados de alfa, beta e gama são 0.41, 0.00 e 0.96, respectivamente. O valor de alfa (0.41) é relativamente baixo, indicando que a estimativa do nível no momento atual é baseada em observações recentes e em algumas observações no passado mais distante. O valor de beta é 0.00, indicando que a estimativa da inclinação b do componente de tendência não é atualizada ao longo da série temporal e, em vez disso, é igual ao seu valor inicial. Isso faz sentido intuitivo, pois o nível muda bastante ao longo da série temporal, mas a inclinação b do componente de tendência permanece aproximadamente a mesma. Em contraste, o valor de gama (0.96) é alto, indicando que a estimativa do componente sazonal no momento atual é baseada apenas em observações muito recentes.

plot(souforecasts)

Vemos no gráfico que o método exponencial de Holt-Winters é muito bem-sucedido em prever os picos sazonais, que ocorrem aproximadamente em novembro de cada ano.

Para fazer previsões para tempos futuros não incluídos na série temporal original procedemos de forma similar ao exemplo anterior.

pred_sou <- forecast(souforecasts,h= 25)
plot(pred_sou)

As previsões são mostradas como uma linha azul, e as áreas sombreadas em cinzamostram intervalos de previsão de 80% e 95%, respectivamente.

acf(pred_sou$residuals, lag.max=20,na.action = na.omit)

Box.test(pred_sou$residuals, lag=20, type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  pred_sou$residuals
## X-squared = 17.53, df = 20, p-value = 0.6183

O correlograma mostra que as autocorrelações para os erros de previsão na amostra não excedem os limites de significância para defasagens 1-20. Além disso, o p-valor para o teste Ljung-Box é 0.6, indicando que há pouca evidência de autocorrelações diferentes de zero nas defasagens 1-20.

Podemos verificar se os erros de previsão têm variância constante ao longo do tempo e são normalmente distribuídos com média zero, fazendo um gráfico de tempo dos erros de previsão e um histograma com curva normal sobreposta

plot.ts(pred_sou$residuals)

plotForecastErrors(pred_sou$residuals)

A partir do gráfico de tempo, parece plausível que os erros de previsão tenham variação constante ao longo do tempo. A partir do histograma de erros de previsão, parece plausível que os erros de previsão sejam normalmente distribuídos com média zero.

Assim, há pouca evidência de autocorrelação nas defasagens 1-20 para os erros de previsão, e os erros de previsão parecem ser normalmente distribuídos com média zero e variância constante ao longo do tempo. Isso sugere que a suavização exponencial de Holt-Winters fornece um modelo preditivo adequado do log de vendas na loja de souvenirs, que provavelmente não pode ser melhorado. Além disso, as suposições sobre as quais os intervalos de previsão foram baseados são provavelmente válidas.

Séries Temporais - Aulas R

Marcio Valk (PPGEST UFRGS)

2022-09-06